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Abstract — The estimation of areal rainfall is an important part of solving various hydrological problems utilizing rainfall - 
runoff and other models. Traditionally used Thiessen polygons (TP) method proved to be inaccurate mainly in mountainous 
areas and in catchments with insufficient number of rainfall gauges. One of the alternative to this method is an inverse 
distance weighting (IDW) method giving better estimates of areal rainfall even on places with rugged orography. However, 
this and similar methods are far less efficient than the simple TP method restraining its use to tasks where computational 
efficiency is not important. In this study a new algorithm accelerating the traditional IDW method is proposed and applied to 
three mountainous catchments situated in the central part of Slovakia. The method is compared with traditional IDW and TP 
methods in terms of both computational efficiency and estimated values. The results showed that while the new method gives 
the same results as the traditional IDW method it is far more efficient when the computational time was in all three 
catchments reduced by more than 96%. 
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I. INTRODUCTION 

In hydrology and other water related disciplines rainfall is one of the most important variables determining the amount of 
water entering the system. Rainfall measurements are essential for solving a wide range of tasks including those in 
meteorology, hydrology, agriculture, climate research or hydropower [7]. Even though the recent years and decades gave rise 
to very sophisticated measurement techniques such as meteorological radars or satellites most of the data is still collected 
with elevated can-type rainfall gauges. These enable to collect precipitation amounts only at single sites. However, rainfall 
data have very strong spatial variability and thus it is not reasonable to presume that a larger area around the rainfall gauge 
will share the same precipitation amounts. In order to account for the spatial variability a network comprising of several rain 
gauges needs to be deployed to correctly represent the rainfall amounts over a certain area. Generally, the denser the network 
is the more accurate is the information about the distribution of rainfall amounts over a certain area. Over the decades a large 
number of models of various complexity have been developed to utilize this information to solve various tasks in hydrology 
and other disciplines. Some of these tasks are of high importance such as predicting floods, managing water resources, 
evaluating crop yields and development or assessing the impact of land use or climate changes on the runoff from a 
catchment. Most of these models are lumped dealing with the catchment as a single unit and thus require estimates of mean 
areal rainfall as model input. The estimation of areal rainfall from a number of rainfall gauges could be performed using a 
number of various techniques. The overview of these techniques and their comparison is given in [1], [5], [8] or [9]. 

One of the most conventional and traditional one is the Thiessen polygon technique [10]. This method is favoured due to its 
simplicity, ease of use and minimum requirements on computational power. However, several authors (see e.g. [6], [8], [2] 
and [3]) comparing the Thiessen polygon technique with other techniques concluded that due to its fundamental principles it 
may produce inaccurate results especially on places with high topographical variation and the limited number of available 
rainfall gauges [9]. One of the possible issues with the Thiessen polygon method is displayed in (Fig. la), where the rainfall 
gauge situated in the upper right corner has no influence on the mean areal rainfall even though it is in the very near vicinity 
of the catchment. 

Another approach in calculating areal rainfall uses inverse distance weighting (IDW) technique belonging to a family of 
distance weighting techniques ([8] and [4]). This technique eliminates the drawback of the Thiessen polygon method by 
dividing the catchment into a grid and interpolating rainfall amounts into each cell of the grid. The technique utilizes data 
from all rainfall gauges while assigning higher weights to stations closer to the interpolated cells. Fig. 2b shows that using 
the IDW technique the upper right station effects the mean areal rainfall over the catchment of interest. One of the 
disadvantages of the method used in calculating areal rainfall is that the precipitation must be interpolated into each cell of 
the grid lying inside the catchment boundaries. Even though the interpolation into particular cells is not computationally 
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demanding in the case of large catchments or dense grids the calculation is repeated numerous times (in case of 500x500 grid 
it is 250 000 times) for each day (or other time step) for which the areal rainfall is calculated. This produces a large number 
of interpolated maps which makes it computationally very demanding (in case of calculating daily areal rainfall amounts for 
10 years this would lead to more than 900 million interpolations). 

When calculating areal rainfall the information about the spatial variability of the precipitation is redundant and not needed. 
By avoiding the interpolation of rainfall amounts into each cell of the grid the computational requirements would decrease 
significantly. In this study an algorithm for accelerating calculation of mean areal rainfall using the IDW method was 
proposed. The new algorithm was compared with traditional IDW method and with the Thiessen polygon method in terms of 
calculation duration and estimated values. 


Fig (A) Fig (B) 



Fig. 1: An example of calculating a real rainfall using two interpolation techniques: a) Thiessen 

POLYGONS AND B) IDW. 

II. Methodology 

Areal rainfall over a certain area A is commonly defined using the following equation: 


P A = — 


,y)dxdy 


( 1 ) 


p,( x >y) 


i*y) 


where denotes the rainfall depth at the point ' for the tth time interval. The total rainfall 1 over the area A is 

unknown, since the rainfall depth is accessible only at finite number of locations n (rainfall gauges). In hydrology the areal 
rainfall is then estimated using so called linear estimators given by 




( 2 ) 


p 1 p 2 p n ^ 

which is a weighted mean of the random variables 1 ’ 1 ’ ' ’ * ’ 1 observed at the rainfall gauges with weights 1 , which sum 
is always 1. Both Thiessen polygons and IDW techniques are linear estimators differing from one another in the values of 

weights w * [5]. 

A. Thiessen polygons method 


In this method the catchment is divided into n zones of influence ' delineating an area which inner points’ closest station is 

area of the 
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catchment A. The weighting coefficients w * are then computed as 

A 


w. = — / = 1, 2, .... ft 
A 


( 3 ) 


The areal rainfall is then estimated from equation (2) which is a sum of weighted rainfall depths from all observed gauges. In 
the case of missing data in one or more gauges the weights had to be recalculated and areal rainfall was estimated from the 
remaining stations. 

B. Inverse distance weighting method 

The IDW interpolation method differs from the Thiessen polygons method in the way the weights are calculated. In the IDW 


method the weights are solely a function of the distances between the point of interest 

(vy,) 


bo>y 0 ) 


and the rainfall gauges 


. This relationship is given by: 


w. =■ 


/( 4 ») 

IfM 


where ^ 0i is the distance between the point of interest and the ith rainfall gauge. Function A^oi) 
given by: 


( 4 ) 

is a distance function 


/k,)=4t 

a 0i 


( 5 ) 


where /? is an appropriate power constant usually from an interval ' ’ ' . Equation (5) indicates that the higher the power 
function is the lower weight is assigned to remote gauges. In this study the power constant was set to be 2. 

When interpolating rainfall depths over large areas or with network of rainfall gauges the information contained in remote 
stations could be neglected since it does not add any relevant additional information to the interpolation procedure (remote 
stations would have negligible weights and thus minor influence on the interpolated value and at places with a sufficient 
amount of observations no additional information is needed). Because of this most of the GIS software contain an additional 
setting limiting the number of stations included to the interpolation either by setting the maximum number of stations or the 
maximum distance from the point of interest (see Fig. 2). 



O 

Fig. 2: An example of selecting stations complying with the maximum distance from the point of interest 

RULE IN THE IDW INTERPOLATION METHOD. 


The estimation of areal rainfall using the traditional IDW technique comprises of several steps. In the first step the rainfall 
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depth is interpolated into all cells of the grid into which the catchment was divided. Since the grid is usually rectangular 
(rainfall depths are interpolated even for cells lying outside the catchment (see Fig. lb)) only cells lying inside the catchment 
are selected. In the last step the mean areal rainfall is calculated by averaging rainfall depths from all cells lying inside the 
catchment (see Fig. 3). In the case of missing data the interpolation is performed from remaining rainfall gauges. 


^ — 


Select next day. 


Interpolate rainfall depths 
into all cells of the grid. 


Select cells inside the 
catchment. 


Calculate areal rainfall. 



Fig. 3: Scheme of the traditional IDW interpolation technique 
C. Accelerated IDW method 

In this section a new algorithm for accelerating the calculation of mean areal rainfall using the IDW method is proposed. This 
method gives a simple ‘shortcut’ to the calculation of areal rainfall using the IDW interpolation technique which would 
otherwise have to be calculated using computationally demanding traditional approach (see the previous section). In order to 
use this technique let’s assume that all rainfall gauges are used in interpolation (no maximum distance and number of stations 
limits). 

P' (jt, y) ( x, y) 

Let v J ' be an interpolated rainfall depth at any arbitrary location with coordinates v 7 ' . This value can be calculated 

P P P W- (i = 1,2, ,n) 

as a linear combination of observed values at stations 1 , 2 , ••• , n and corresponding weights 1 v } 

Rainfall depth at any location can be calculated by 


p \x,y) = Y j w i P i 

( 6 ) 

p r ( x, y ) 

Consider, the problem of finding the average within a catchment A. This average will be found by predicting ' 
separately for each grid cell within the catchment and averaging those predictions. Let the cell centres have coordinates 

(( w i’ v i )’ • • v m )) w j iere m j s th e num ber of cells in the region. The average ^ is defined as 


m 

P'(A) = ^ 

m 

By combining equations (6) and (7) and their subsequent modification we can get: 


(7) 


p\A = ±p;~ = Z p r w 'i 


m 


( 8 ) 


where w ‘ is only an average weight of a particular station over the whole catchment A. Let us recall that the individual 
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rainfall gauges contain a time series of observed rainfall depths. Let us denote rainfall depth at station i and time t u . Since 
the location of the gauges and the size of the grid cells is constant even the weights W ‘ doesn’t change in time. This implies 

w' P 

that the areal rainfall in time t can be calculated as a linear combination of weights 1 and observed rainfall depths u . 

P w' 

Since the rainfall depths at stations u are known only the n weights 1 have to be estimated. This can be done by setting 

p 

value at station u at 1 and in all other stations at 0. After substituting these values to equation (8) we get 

P' ( A) = Wj • 0 + . . . + w'- • 1 + w^ +1 • 0 + . . . + w’ n • 0 

p'(a)= w ; (9) 


Equation (9) implies that the weights w * can be calculated as average values of synthetic rainfall depths (1 for i, 0 for others) 
over the catchment A. The calculation of the mean areal rainfall could be then reduced into the following steps: 

Step 1: Set ^ and ^ kt ~ ^ for ^ ^ * . 

• Interpolate these values into each cell of the grid. 

• Calculate mean of the interpolated values for each cell lying inside the catchment to get W/ . 


Step 2: Repeat step 1 for each station to get all weights 


w. 


P w 

Step 3: Calculate areal rainfall as a linear combination of observed rainfall depths 11 and weights * 
Step 4: Repeat step 3 for each day with observations. 

^ Start 
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missing stations. 
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database. 




Fig. 4: Scheme of the accelerated IDW interpolation technique. 
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In the case of missing data in one or more stations it is necessary to recalculate the weights with the remaining stations 
(stations with missing data are not included into the interpolation). Since the same combination of missing stations could 
occur anywhere in the observed data (the same weights would have to be calculated several times) each combination of 
weights is stored in the weight database. If some of the stations at day t are missing the algorithm searches the database for 
the right combination of missing stations. If such combination is not found the new weights are calculated and stored in the 
database. The scheme of the algorithm is displayed in (Fig. 4). 


III. Data 


The estimation of mean daily areal rainfall was tested on three mountainous catchments situated in the central part of 
Slovakia. The catchments are: 1) Hron at Banska Bystrica, 2) Hron at Brezno and 3) Vah at Liptovky Mikulas (for the 
position of the catchments see Fig. 5). In all the catchments daily rainfall data were collected for the period between 1.1.1961 
and 31.12.2010. 



Fig. 5: Situation of catchments of interest in the centre of Slovakia. 

Catchment Hron with the outlet at Banska Bystrica (Hron@BB) has an area of 1768.2 km 2 . The catchment has a shape of a 
narrow valley with two high mountains situated in the northern and southern part of the catchment. The average altitude of 
the catchment is 845 m ASL with minimum and maximum point at 350 and 2043 mASL respectively. The areal rainfall was 
estimated from 24 rainfall gauges situated inside the catchment and in its near vicinity. 

Catchment Hron with the outlet at Brezno (Hron@Brezno) is a nested catchment of Hron@BB. It has an area of 578.6 km 2 . 
Its mean altitude is 916 m ASL with minimum and maximum elevation of 489 and 1946 m ASL respectively. Rainfall depths 
were collected from 1 1 rainfall gauges situated mostly in the lower parts of the catchment. 

Catchment Vah with the outlet at Liptovsky Mikulas (Vah@LM) is a neighbouring catchment to the Hron@BB with which it 
shares the southern border. Its area is 1 100.6 km 2 with the mean altitude of 1090 m ASL (min. 568 m ASL and max. 2494 m 
ASL). Rainfall depths were collected from 28 rainfall gauges with some of the gauges situated in the highest parts of the 
catchment. 

In all three catchments the rainfall datasets contain a certain percentage of missing data. The percentage of days with at least 
one missing value in any station and the percentage of missing values from all data is described in (Table 1). 

Table 1 

Basic characteristics of selected catchments. The last two columns describe the percentage of days 

WITH AT LEAST ONE MISSING DATA IN ANY STATION AND PERCENTAGE OF ALL MISSING VALUES. 


Catchment 

Area 

[km 2 ] 

Mean alt. 
[m ASL] 

Number of 
gauges 

Period 

Days with 
missing data 

Missing 

values 

Hron@BB 

1768.2 

845 

24 

1961 -2010 

83% 

9.1% 

Hron @ Brezno 

578.6 

916 

11 

62% 

9.6% 

Vah@LM 

1100.6 

1090 

28 

77% 

29% 
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IV. Results and Discussion 

In order to assess the performance of the new algorithm (IDW-a) the calculation of mean areal rainfall was compared with 
the traditional IDW method (IDW-t) and with the Thiessen polygons method (TP). The methods were compared both in 
terms of the estimated areal rainfall and their computational requirements. 

A. Comparison of estimated areal rainfall 

Since the main objective of the study was to propose a method which would accelerate the calculation of areal rainfall using 
the tradional IDW method it was crutial to demonstrate that the new method gives identical results. The method of Thiessen 
polygons was also included into the comparison to prove that it gives different estimates of areal rainfall. The individual 
methods were compared to the traditional IDW method (IDW-t) which was used as a reference method. The comparison 
showed that the difference between the estimates of areal rainfall calculated using the IDW-a and the IDW-t methods are 
negligible and thus could be considered as identical (see Table 2). When comparing the TP to the reference IDW-t method it 
is possible to observe a small difference between these two methods. Even though the maximum absolute difference between 
estimeted areal rainfall is almost 4 mm and in the case of Vah@LM more than 7 mm the mean absolute error is in all 
catchments lower that 0.2 mm. This means that even though the TP method doesn’t usually give very accurate estimates in 
this particular case the estimates are comparable to those calculated by the IDW methods. 


Table 2 

The comparison of selected statistics between TP and IDW-t methods and IDW-a and IDW-t 

METHODS 


Catchment 

Hron@BB 

Hron@Brezno 

Vah@LM 


TP 

IDW-a 

TP 

IDW-a 

TP 

IDW-a 

Max abs. err. 

3.9 

<0.001 

3.8 

<0.001 

7.3 

<0.001 

MAE 

0.111 

<0.001 

0.123 

<0.001 

0.187 

<0.001 

RMAE 

0.034 

<0.001 

0.34 

<0.001 

0.059 

<0.001 


B. Comparison of algorithm efficiency 

The computational efficiency of the methods in estimating mean areal rainfall was tested at two datasets for each catchment. 
The first dataset comprised of 50 years of daily rainfall depths observed at individual rainfall gauges during the period 
between 1961 and 2010. This dataset was not complete and contained a certain percentage of missing data in particular 
stations (see Table 1). The computational efficiency of the methods was evaluated based on the time required for the 
estimation of areal rainfall in all days of the dataset. The results of the analysis are shown in Fig. 6 and Table 3. In all three 
catchments the slowest method is IDW-t method followed by IDW-a and TP method. In the Hron@BB catchment the time 
required for the estimation of areal rainfall using the traditional IDW-t method is over 5 hours and 43 minutes. The identical 
values (see Table 2) of areal rainfall were estimated using the proposed accelerated IDW-a method in less than 15 minutes 
representing a 96% reduction in computational time. A similar pattern could be observed in all other catchments when the 
computational time was reduced by 97% in Vah@LM and 99% in Hron@Brezno. As expected the highest computational 
efficiency was observed in the case of the TP method ranging only around 1 minute. Since the accuracy of the TP method is 
questionable especially in mountainous catchments [9] and the objective of the study was only to accelerate the estimation of 
areal rainfall using the IDW method the discussion was foccused on the comparison of IDW-t and IDW-a methods 



Fig. 6: Comparison of time needed for estimation of areal rainfall for the whole dataset with missing 

DATA. 
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Table 3 

Time needed for calculation of mean daily areal rainfall using the Thiessen polygon (TP), traditional 
IDW (IDW-t) and accelerated IDW (IDW-a) method. Calculation for the period between 1961 and 2010 

WITH MISSING DATA. 


Catchment 

Elapsed time 


TP 

IDW-t 

IDW-a 

Hron@BB 

70s 

5h 43m 11s 

14m 22s 

Hron@Brezno 

44s 

4h 26m 56s 

2m 20s 

Vah@LM 

59s 

5h 12m 17s 

9m 43 s 


One of the task with the highest computational requirements in the IDW-a method is the calculation of weights. This process 
requires to interpolate synthetic rainfall depths into each cell of the grid into which the catchment is divided. This process is 
repeated for each rainfall gauge and for each observed combination of missing stations. In the case of a large number of 
missing data in an analyzed dataset the new IDW-a method could be less efficient and comparable with traditional IDW-t 
method. In order to test the algorithm in ideal conditions in which the dataset would be complete (no missing data) a new test 
was conducted. In this test all the three methods were used to estimate areal rainfall over the three catchments for a year with 
no missing data. The time requirements of individual methods and catchments are displayed in Fig. 7 and listed in Table 4. 
The results are comparable with those from the previous test when the whole dataset with missing data was used. In all three 
catchments the time needed for the estimation of areal rainfall using the TP method was around 2 seconds (see Table 4) 
which is less than 1% of the time needed when using the IDW-t method (reference method). The time requirements for the 
IDW-a method constituted only for 3 to 7% of the IDW-t method. 



Fig. 7: Comparison of time needed for estimation of areal rainfall for a period of 1 year without missing 

DATA. 


Table 4 

Time needed for calculation of mean daily areal rainfall using the Thiessen polygon (TP), traditional 
IDW (IDW-t) and accelerated IDW (IDW-a) method. Calculation for 1 year without missing data. 


Catchment 

Elapsed time 


TP 

IDW-t 

IDW-a 

Hron@BB 

2s 

7min 12s 

29s 

Hron@Brezno 

1.8s 

5min 25 s 

10s 

Vah@LM 

2s 

6min 48 s 

30s 


V. Conclusion 

The process of estimation areal precipitation is of high importance and is an integral part of solving various tasks in 
hydrology and other disciplines utilizing mainly conceptual rainfall -runoff models which are not distributed but work as 
lumped models. The number of methods for estimation of areal rainfall differ mainly in their accuracy but also in 
computational requirements. One of the methods which proved to be more accurate than the simple method utilizing the 
Thiessen polygons is an inverse distance interpolation method [8]. However, this method is computationally far less effective 
than the TP method. The aim of this study was to propose an algorithm which would significantly reduce the time needed for 
estimation of areal rainfall using the traditional IDW-t method. The new algorithm was tested not only in the terms of 
computational efficiency but also in the values areal rainfall it estimated. The algorithm was compared with the TP and IDW - 
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t methods on data from three mountainous catchments situated in the central part of Slovakia. Two tests on two different 
datasets were performed for each catchment: 1) the whole dataset with missing data and 2) one year with no missing data. 

The results showed that when comparing the values of estimated areal rainfall the IDW-a and IDW-t methods give practically 
identical results with mean absolute error in all three catchments lower than 0.001 (Table 2). As expected the estimates of the 
TP method when compared to the IDW-t method were slightly different with the maximum absolute error ranging from 3.8 
mm (Hron@Brezno) to 7.3 mm (Vah@LM). When comparing the computational efficiency of the methods the new IDW-a 
method has significantly reduced the time needed for the estimation of areal rainfall using the traditional IDW -t method. In 
all three catchments the calculation was reduced by more than 96% in the case of analyzing the whole dataset (Table 3) and 
by more than 92% in the case of analyzing only 1 year (Table 4). In the case of the Hron@Brezno catchment the 
computational time was reduced by more than 99% when with the traditional IDW-t method it took 4h 26m 56s and with the 
new IDW-a method only 2m 20s. Similar results were obtained in the remaining catchments. 

The significant reduction of computational time in the new IDW-a method could be used in solving many practical tasks. 
One of them is preparing areal rainfall data for zonal conceptual rainfall -runoff models where for each zone (most often 
based on altitude) a separate dataset of rainfall depths is prepared. The computational efficiency of the new method enables 
to divide the catchment into various number of zones and select the best division based on the selected criteria (e.g. 
comparison of observed and simulated flows, simulation of accumulation and melting of snow cover). 

The process of accelerating the IDW method described in this study could be also applied to other linear estimators such as 
Kriging. 
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